This article demonstrates how to work with coastal catchments and small coastal basins in an hyRefactor workflow.
First, we use nhdplusTools to download and plot some data along the California coast.
library(hyRefactor) #> USGS Support Package: https://owi.usgs.gov/R/packages.html#support temp_gpkg <- tempfile(fileext = ".gpkg") bbox <- sf::st_bbox(c(xmin = -124.4, ymin = 39.8, xmax = -123.7, ymax = 40.3), crs = sf::st_crs(4326)) nhd_data <- nhdplusTools::plot_nhdplus(bbox = bbox, gpkg = temp_gpkg, overwrite = TRUE) #> Zoom: 10 #> Map tiles by Carto, under CC BY 3.0. Data by OpenStreetMap, under ODbL. #> Audotdetect projection: assuming Google Mercator (epsg 3857)
Now that we have some data to work with, we need to subset it so we have a few complete basins to work with.
Note that both basins that terminate at the coast AND the coastal flowlines are included in out subset below.
First we use the nhdplusTools::navigate_nldi for each terminal flowline in our subset to find complete basins then we use the nhdplusTools::subset_nhdplus to get the data we need.
terminals <- dplyr::filter(nhd_data$flowline, TerminalFl == 1) coastal <- dplyr::filter(nhd_data$flowline, FTYPE == "Coastline") flowline <- dplyr::bind_rows( lapply(terminals$COMID, function(x) { nhdplusTools::navigate_nldi(list(featureSource = "comid", featureID = x), mode = "UT") })) sub <- lapply(nhdplusTools::subset_nhdplus(c(flowline$nhdplus_comid, coastal$COMID), nhdplus_data = "download"), nhdplusTools::align_nhdplus_names) #> All intersections performed in latitude/longitude. #> Reading NHDFlowline_Network #> Writing NHDFlowline_Network #> Reading CatchmentSP #> Writing CatchmentSP mapview::mapview(sub$NHDFlowline_Network) + mapview::mapview(sub$CatchmentSP, color = "darkgrey", col.regions = "tan")
fdr_dir = file.path(tempdir(check = TRUE), "fdr_fac") dir.create(fdr_dir, showWarnings = FALSE) if(!dir.exists(file.path(fdr_dir, "NHDPlusCA/"))) { hyRefactor::download_fdr_fac(fdr_dir, regions = "18") } fdr_path <- file.path(fdr_dir, "NHDPlusCA/NHDPlus18/NHDPlusFdrFac18c/fdr/") fac_path <- file.path(fdr_dir, "NHDPlusCA/NHDPlus18/NHDPlusFdrFac18c/fac/")
fdr <- raster::raster(fdr_path) fac <- raster::raster(fac_path) crs <- raster::crs(fac) fdr <- raster::crop(fdr, sf::as_Spatial(sf::st_transform(sub$CatchmentSP, crs))) fac <- raster::crop(fac, sf::as_Spatial(sf::st_transform(sub$CatchmentSP, crs))) mapview::mapview(fdr) #> Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; #> the supplied Raster* has 3148995 #> ... decreasing Raster* resolution to 5e+05 pixels #> to view full resolution set 'maxpixels = 3148995 ' #> Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : Using #> PROJ not WKT2 strings #> Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : Using #> PROJ not WKT2 strings #> Warning in rgdal::rawTransform(projto_int, projfrom, nrow(xy), xy[, 1], : Using #> PROJ not WKT2 strings #> Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : Using #> PROJ not WKT2 strings #> Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : Using #> PROJ not WKT2 strings #> Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : Using #> PROJ not WKT2 strings #> Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : Using #> PROJ not WKT2 strings #> Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : Using #> PROJ not WKT2 strings #> Warning in rgdal::rawTransform(projto_int, projfrom, nrow(xy), xy[, 1], : Using #> PROJ not WKT2 strings
net <- sub$NHDFlowline_Network cats <- sub$CatchmentSP coastal <- net[net$FTYPE == "Coastline", ] coastal_cats <- cats[cats$FEATUREID %in% coastal$COMID, ] net <- dplyr::filter(net, !COMID %in% coastal$COMID) cats <- cats[cats$FEATUREID %in% net$COMID, ] nhd_outlets <- dplyr::filter(net, TerminalFl == 1) sf::st_geometry(nhd_outlets) <- sf::st_geometry(nhdplusTools::get_node(nhd_outlets, "end")) mapview::mapview(net) + mapview::mapview(cats, color = "darkgrey", col.regions = "tan") + mapview::mapview(coastal, color = "lightblue") + mapview::mapview(coastal_cats, color = "darkgrey", col.regions = "brown") + mapview::mapview(nhd_outlets, color = "black")
Now we will identify coastal basins that we want to lump into coastal catchments as part of the refactor workflow.
min_da_km <- 10 little_terminal <- dplyr::filter(net, TerminalPa %in% dplyr::filter(nhd_outlets, TotDASqKM <= min_da_km & TerminalFl == 1)$TerminalPa) outlets <- dplyr::select(nhd_outlets, COMID) %>% dplyr::mutate(type = "terminal") %>% dplyr::filter(COMID %in% cats$FEATUREID) %>% dplyr::mutate(keep = ifelse(COMID %in% little_terminal$COMID, "temporary", "keep")) mapview::mapview(outlets, zcol = "keep") + mapview::mapview(net)
We can now run refactor_nhdplus and reconcile_catchment_divides.
Note that the exclude_cats parameter is set to all outlet flowlines in small basins that were identified above.
Also note that the net variable no longer contains coastal flowlines.
tf <- file.path(tempfile(fileext = "tf.gpkg")) tr <- file.path(tempfile(fileext = "tr.gpkg")) refactor_nhdplus(nhdplus_flines = net, split_flines_meters = 100000, split_flines_cores = 1, collapse_flines_meters = 2000, collapse_flines_main_meters = 2000, out_refactored = tf, out_reconciled = tr, three_pass = TRUE, purge_non_dendritic = FALSE, exclude_cats = unique(c(outlets$COMID, little_terminal$COMID)), warn = FALSE) refactored <- sf::st_transform(sf::read_sf(tf), crs) reconciled <- sf::st_transform(sf::read_sf(tr), crs) cats <- sf::st_transform(cats, crs) sf::st_precision(cats) <- 10 divides <- reconcile_catchment_divides(catchment = cats, fline_ref = refactored, fline_rec = reconciled, fdr = fdr, fac = fac, para = 1) mapview::mapview(reconciled) + mapview::mapview(divides)
keep_outlets <- dplyr::filter(outlets, keep == "keep") %>% dplyr::select(COMID, type) mapped_outlets <- map_outlet_ids(keep_outlets, reconciled) %>% dplyr::filter(COMID %in% keep_outlets$COMID) zero_order <- list(basin = little_terminal$COMID, zero = coastal$COMID) agg_cats <- aggregate_catchments(flowpath = reconciled, divide = divides, outlets = dplyr::select(mapped_outlets, ID, type), zero_order = zero_order, coastal_cats = sf::st_transform(coastal_cats, sf::st_crs(divides)), da_thresh = 1, only_larger = TRUE) #> removing 2 headwater/diversion flowlines without catchments. mapview::mapview(agg_cats)